| age | lx | dx | mx | qx |
|---|---|---|---|---|
| 0 | 100000 | 553.22 | 0.00556 | 0.00553 |
| 1 | 99447 | 38.18 | 0.00038 | 0.00038 |
| 2 | 99409 | 23.16 | 0.00023 | 0.00023 |
| age | lx | dx | mx | qx |
|---|---|---|---|---|
| 0 | 100000 | 553.22 | 0.00556 | 0.00553 |
| 1 | 99447 | 38.18 | 0.00038 | 0.00038 |
| 2 | 99409 | 23.16 | 0.00023 | 0.00023 |
lx is the total number living in a hypothetical cohort (of 100,000).dx is the number of deaths in the age interval.| age | lx | dx | mx | qx |
|---|---|---|---|---|
| 0 | 100000 | 553.22 | 0.00556 | 0.00553 |
| 1 | 99447 | 38.18 | 0.00038 | 0.00038 |
| 2 | 99409 | 23.16 | 0.00023 | 0.00023 |
mx is the death rate within the interval.qx is the probability of death within an interval.| age | lx | dx | mx | qx |
|---|---|---|---|---|
| 0 | 100000 | 553.22 | 0.00556 | 0.00553 |
| 1 | 99447 | 38.18 | 0.00038 | 0.00038 |
| 2 | 99409 | 23.16 | 0.00023 | 0.00023 |
| age | lx | dx | mx | qx |
|---|---|---|---|---|
| 0 | 100000 | 553.22 | 0.00556 | 0.00553 |
| 1 | 99447 | 38.18 | 0.00038 | 0.00038 |
| 2 | 99409 | 23.16 | 0.00023 | 0.00023 |
The R packages MortalityLaws and demography facilitate life table data aquisition and modeling.
Common sources:
Highlighted column: Remaining life expectancy at exact age x.
| age | lx | dx | ex | qx | mx |
|---|---|---|---|---|---|
| 0 | 100000 | 553.221 | 78.963 | 0.00553 | 0.00556 |
| 1 | 99447 | 38.180 | 78.402 | 0.00038 | 0.00038 |
| 2 | 99409 | 23.160 | 77.432 | 0.00023 | 0.00023 |
| 3 | 99385 | 17.689 | 76.450 | 0.00018 | 0.00018 |
qx) values.params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)lifetable “strategy” will draw on the HMD life table data processed earlier.We’ll next construct a transition probability matrix using the age-specific probability of death field (qx) from the life table data.
To do so, we need to define a function that, for a given age and cycle length, calculates the probability of death and places this within a transition matrix.
Life Table:
sim_cohort <- function(params) {
params$t_names %>% map(~({
tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
res <- do.call(rbind,lapply(params$mP, function(tp) {
tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
}))
res <- rbind(c(params$n_cohort,0),res)
dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
res
})) %>%
set_names(params$t_names)
}| age | qx | ex |
|---|---|---|
| 0 | 0.00553 | 78.963 |
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)| Age | Life Expectancy (Life Table data) | Life Expectancy (Markov-Life Table) |
|---|---|---|
| 0 | 78.963 | 78.924 |
MortalityLaws package has a number of mortality models we can draw from:| Name | Model | Code |
|---|---|---|
| Infant mortality | ||
| Opperman | \(\mu[x] = A/\sqrt{x} - B + C\sqrt{x}\) | opperman |
| Weibull | \(\mu[x] = 1/\sigma * (x/M)^{(M/\sigma - 1)}\) | weibull |
| Accident hump | ||
| Inverse-Gompertz | \(\mu[x] = [1-\exp((M-x)/\sigma)]/[\exp((M-x)/\sigma)-1]\) | invgompertz |
| Inverse-Weibull | \(\mu[x] = 1/\sigma * (x/M)^{[-M/\sigma - 1]} / [\exp((x/M)^{(-M/\sigma)}) - 1]\) | invweibull |
| Name | Model | Code |
|---|---|---|
| Adult mortality | ||
| Gompertz | \(\mu[x] = A \exp[Bx]\) | gompertz |
| Gompertz | \(\mu[x] = 1/\sigma \exp[(x-M)/\sigma]\) | gompertz0 |
| Makeham | \(\mu[x] = A \exp[Bx] + C\) | makeham |
| Makeham | \(\mu[x] = 1/\sigma \exp[(x-M)/\sigma)] + C\) | makeham0 |
| Perks | \(\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]\) | perks |
| Strehler-Mildvan | \(\mu[x] = K \exp[-V_0 (1-Bx)/D]\) | strehler_mildvan |
| Adult and/or old-age mortality | ||
| Van der Maen | \(\mu[x] = A+Bx+Cx^2 + I/[N-x]\) | vandermaen |
| Beard | \(\mu[x] = \exp(Bx)/[1+KA \exp(Bx)]\) | beard |
| Beard-Makeham | \(\mu[x] = A \exp(Bx)/[1 + KA \exp(Bx)]\) | beard_makeham |
| Gamma-Gompertz | \(\mu[x] = A \exp(Bx)/(1+AG/B[\exp(Bx)-1])\) | ggompertz |
| Name | Model | Code |
|---|---|---|
| Old-age mortality | ||
| Van der Maen | \(\mu[x] = A + Bx + I/[N-x]\) | vandermaen2 |
| Quadratic | $= A + Bx + Cx^2% | quadratic |
| Kannisto | \(\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]\) | kannisto |
| Kannisto-Makeham | \(\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]+C\) | kannisto_makeham |
| Name | Model | Code |
|---|---|---|
| Full age range | ||
| Thiele | \(\mu[a] = A e^{-Bx} + C e^{-\frac{1}{2} D(x-E)^2}+F e^{Gx}\) | thiele |
| Wittstein | \(q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}\) | wittstein |
| Heligman-Pollard | \(q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + G H^x\) | HP |
| Heligman-Pollard | \(q[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + GH^x / [1 + GH^x]\) | HP2 |
| Rogers-Planck | \(q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}\) | rogersplanck |
| Martinelle | \(\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}\) | martinelle |
| Carriere | \(l[x] = P1 l[x](weibull) + P2 l[x](invweibull) + P3 l[x](gompertz)\) | carriere1 |
| Kostaki | \(q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x\) | kostaki |
Generally speaking, we need three inputs:
age: Ages for lifetabledx: The number of deaths between exact ages x and x+1.lx: Number of survivors to exact age x. A B
0.00010772 0.07534542
# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 ) A
0.0021938
# Can simply use the supplied function to calcualte
MortalityLaws::gompertz(x = 40, par = gom_fit$coefficients)$hx
[1] 0.0021938
$par
A B
0.00010772 0.07534542
$Sx
[1] 0.97269
ages <- lt$age[lt$age>=40 & lt$age<100]
deaths <- lt$dx[lt$age>=40 & lt$age<100]
exposure <- lt$lx[lt$age>=40 & lt$age<100]
gom_fit40 <- MortalityLaw(
x = ages, # vector with ages
Dx = deaths, # vector with death counts
Ex = exposure, # vector containing exposures
law = "gompertz",
opt.method = "LF2")While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP in lieu of gompertz!
While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP in lieu of gompertz!
Mortality rate for a 40 year old:
$hx
[1] 0.0020409
$par
A B C D E F_
0.000386408 0.021271772 0.107422357 0.001025285 2.871760263 33.201213579
G H
0.000022277 1.102611045
[1] 70.467
| age | mx | Sx |
|---|---|---|
| 70 | 0.01827 | 0.76590 |
| 71 | 0.01995 | 0.75077 |
| 72 | 0.02162 | 0.73472 |
y for a given value of x.approxfun():
x is quantile in CDF for survival, i.e., 1 - Survival1y is the age.qHP <-
approxfun(y = 0:115,
x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(0.25)[1] 71.069
Survival in life table:
| age | mx | Sx |
|---|---|---|
| 70 | 0.01827 | 0.76590 |
| 71 | 0.01995 | 0.75077 |
| 72 | 0.02162 | 0.73472 |
params =
list(
t_names = c("lifetable", # Strategy names
"heligman-pollard",
"gompertz"),
n_treatments = 3, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 40, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt, # Processed HMD life table data (2019, USA)
gom_fit = gom_fit, # Fitted Gompertz model object
hp_fit = hp_fit # Fitted Heligman-Pollard model object
)fn_mPt <- function(t, params) {
with(params, {
h = cycle # Time step
lapply(t, function(tt){
# Life table method
current_age_lt = pmin(initial_age + (tt)*h - 1,max(lt$age))
p_death <- lt[lt$age==current_age_lt,"qx"]
# parametric mortality methods
current_age = initial_age + (tt)*h - 1
r_death_gompertz <- gompertz(current_age,params$gom_fit$coefficients)$hx
p_death_gompertz <- 1 - exp(-r_death_gompertz)
r_death_hp <- HP(current_age, params$hp_fit$coefficients)$hx
p_death_hp <- 1 - exp(-r_death_hp)
array(data = c(1 - p_death, 0,
p_death,1,
1 - p_death_hp, 0,
p_death_hp,1,
1 - p_death_gompertz, 0,
p_death_gompertz,1),
dim = c(n_states, n_states, n_treatments),
dimnames = list(from = s_names,
to = s_names,
t_names))
})
})
}# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = lapply(trace, function(tr) (tr * (1 / params$n_cohort)) %*% payoff_life_exp)
# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = lapply(life_exp_cycle,cycle_adj)Constructing the cause-deleted life table requires two necessary inputs:
We extract the percentage of overall deaths that are from CVD by age bin:
ihme_cvd <-
tibble::tribble(
~age_name, ~val,
1, 0.038771524,
5, 0.038546046,
10, 0.044403585,
15, 0.033781126,
20, 0.035856165,
25, 0.053077797,
30, 0.086001439,
35, 0.130326551,
40, 0.184310334,
45, 0.21839762,
50, 0.243705394,
55, 0.256334637,
60, 0.26828001,
65, 0.272698709,
70, 0.28529754,
75, 0.310642009,
0, 0.016750489,
80, 0.353518012,
85, 0.399856716,
90, 0.447817792,
95, 0.495305502
) %>%
mutate(age_ihme = cut(age_name,unique(c(0,1,seq(0,95,5),105)),right=FALSE)) %>%
select(age_ihme, pct_cvd = val) Merge these percentages into the underlying life table data and use them to calculate the total number of CVD deaths by age:
| age | age_ihme | pct_cvd | dx | dx_i |
|---|---|---|---|---|
| 0 | [0,1) | 0.01675 | 553.221 | 9 |
| 10 | [10,15) | 0.04440 | 11.518 | 1 |
| 25 | [25,30) | 0.05308 | 103.706 | 6 |
| 50 | [50,55) | 0.24371 | 365.721 | 89 |
| 75 | [75,80) | 0.31064 | 1997.766 | 621 |
| 98 | [95,105) | 0.49531 | 1262.593 | 625 |
We next calculate the cause-deleted probability of death (qd) and death rate (md), as well as the cause-specific probability of death (qi) and death rate (mi):
qi): \(q \cdot dx_i / dx\)qi): \(q \cdot dx_i / dx\)| age | dx | dxd | dx_i | q | qi | md | mi |
|---|---|---|---|---|---|---|---|
| 0 | 553.221 | 544.221 | 9 | 0.00553 | 0.00009 | 0.00546 | 0.00009 |
| 1 | 38.180 | 37.180 | 1 | 0.00038 | 0.00001 | 0.00037 | 0.00001 |
| 2 | 23.160 | 22.160 | 1 | 0.00023 | 0.00001 | 0.00022 | 0.00001 |
| 3 | 17.689 | 16.689 | 1 | 0.00018 | 0.00001 | 0.00017 | 0.00001 |
| 4 | 14.209 | 13.209 | 1 | 0.00014 | 0.00001 | 0.00013 | 0.00001 |
| 5 | 13.213 | 12.213 | 1 | 0.00013 | 0.00001 | 0.00012 | 0.00001 |
ages_ <- lt_$age[lt_$age<100 & lt_$age>=25]
deaths_ <- lt_$dxd[lt_$age<100 & lt_$age>=25]
exposure_ <- lt_$lx[lt_$age<100 & lt_$age>=25]
mort_fit_CVDdeleted <- MortalityLaw(
x = ages_,
Dx = deaths_, # vector with death counts
Ex = exposure_, # vector containing exposures
law = "HP2",
opt.method = "LF2") Healthy CVD CVDDeath Dead
0 100000 0.0 0.0000 0.00
1 90391 9506.2 0.2945 103.00
2 81702 18088.0 1.1404 208.58
3 73847 25833.6 2.4853 316.88
4 66745 32822.6 4.2817 428.11
5 60324 39126.6 6.8544 542.45
6 54519 44809.4 11.5751 660.14
7 49270 49930.7 17.3944 781.42
8 44526 54544.1 23.8195 906.57
9 40236 58697.0 31.3741 1035.87